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Correlation of macroscopic instability and Lyapunov times 
in the general three-body problem 
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ABSTRACT 

We conducted extensive numerical experiments of equal mass three-body systems until they 
became disrupted. The system lifetimes, as a bound triple, and the Lyapunov times show a 
correlation similar to what has been earlier obtained for small bodies in the Solar System. 
Numerical integrations of several sets of differently randomised initial conditions produced 
the same relationship of the instability time and Lyapunov time. Marginal probability densities 
of the various times in the three-body experiments are also discussed. Our high accuracy 
numerical method for three-body orbit computations and Lyapunov time determinations is 
concisely described. 
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1 INTRODUCTION 

iLecar. Franklin. & Murisor] i 19921) first suggested the existence of 
a relationship between the local Lyapunov time and a survival time 
of the orbit of an asteroid in the Solar System. They gave numerical 
evidence of a powerlaw relation 
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(1) 



where td is the survival time, t e the local Lyapunov time [which 
is defined as the inverse of the mean time derivative of the loga- 
rithm of the variational equation solution during the interplay pe- 
riod of the system] and to some normalisation time. For the pow- 
erlaw index they estimated the value n ~ 1.8. The relation is 
by no means an exact one but rather the equation fits the centre 
of a scatter diagram of the logarithms of the times. The relation 
have been studied by several authors but mainly for asteroid mo- 
tions and for s ome artificial mappings. Additiona l dis cussion can 
be fou n d e.g. in Murison. Lecar. & Franklin! dl994h and lLecar et al.1 
feOOlh . lMorbidelli & Froeschlel dl996h discovered two separate re- 
gions: the Nekhoroshev region in which the survival times are ex- 
ponentially long, and resonance overlap region in which an expres- 
sion like Eq.lQJ more accurately apply. Howev er, they were unable 
to pro vide any theoretical explanation. Later Murr ay & HolmarJ 
i 19971) obtained evidence against the generality of the 'law' iQJ. 

In this paper we present evidence for a relation of the power 
law type in the general equal mass three-body problem. 



2 PROBLEM AND DEFINITIONS 

The problem studied in this paper is the relationship between life- 
times and Lyapunov times of a three-body system. 

The lifetime is defined as the duration of the numerical experi- 
ment to the point when the system ejects one particle. The criterion 
used was: the single particle is moving away in a hyperbolic relative 
orbit at a distance > 50 times the semi-major-axis of the surviving 
binary. Then, using the orbital elements, we compute the (formal) 
time of pericentre passage for this relative hyperbola. This time, 
symbol td, is then considered to be the moment of disruption. 

In the strict mathematical sense the Lyapunov time does not 
exist for unstable triples, because after the systems disrupt they 
become essentially quasi-periodic. We can still use the variational 
equations solutions to measure the sensitivity of the system with 
respect to variations of initial conditions. If the system behaves 
chaotically, then Sr(t), the solution of the variational equations, 
has the order of magnitude |<5r(£)| ~ exp(i/t e ). Evaluating this 
at the time t = td and solving for t e gives (formally) a Lyapunov 
time as 

t a = t d /\n\5r(t d )\. (2) 
Writing the above in another way, one notes that the quantity 

Z = ln\5r(t d )\ =t d /t e , (3) 

i.e. the logarithm of the variation at the moment of disruption, is 
actually the lifetime in units of the Lyapunov time. 
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Figure 1. Labelling of vectors in the three-body regularization. 



3 NUMERICAL METHOD 

Our n umerical method is based on the logarithmic Hamilto - 
nian jMikkola & Tanikaw3 Jl999al lbh: IPreto & Tremaind dl999h ; 
Mikkola (2006)) which, with the leapfrog algorithm, gives regu- 
lar results. Instead of centre-of-mass coordinates we use the three 
interparticle vectors (see Figure Q} 



Xi=r 3 -r 2 ; X 2 — ri - r 3 ; X 3 =r 2 -ri. 



(4) 



as new coordinates. Let the corresponding velocities be V k = X k, 
then the centre-of-mass kinetic and potential energies can be writ- 
ten 



i<3 



miinj 



(5) 



where M = m k is the total mass and kij = 6 
The equations of motion are 

X k = V k ; V k = -M^ +mk J2j^ 



(6) 



and after the application of the logarithmic Hamiltonian modifica- 
tion they read 

t' = l/(T + B); X' k = X k /(T + B); V' k = V k /U. (7) 

where the constant B = U — T is the negative of energy. The 
leapfrog algorithm is now possible since the right hand sides in 
these equations do not depend on the left hand side variable. 

To obtain the solutions to the variational equations 
8X,5t, SV we first wrote the code to compute the motion and 
then differentiated that line by line (i.e. adding the (analytical) code 
lines for the differentials), thus obtaining precise differentials of the 
numerical algorithm used. The results are then in the time trans- 
formed system in which time also has a variation. The transforma- 
tion to physical system is simple: If we denote the physical system 
variations with A, we can write for any quantity / the equation 



Af = 5f-Stf, 



(8) 



which gives the desired result. 

Leapfrog is used only as a basic method and the results are 
improved to high precisio n by means of the extrapolation method 
dBulirsch and Stoerl < fl966h ). 

The usage of the relative vectors, instead of some inertial 
coordinates, is advantageous in attempting to avoid large round- 



off effects. One could also integrate only two of the triangle 
sides, obtaining the remaining one from the conditions J^ fc X k = 
0; Y^ k = 0. This would hardly reduce the computational ef- 
fort required by the method. Instead one may, occasionally, com- 
pute the longest side, and the corresponding velocity, from the 
above triangle conditions. Note that the sums of the sides are not 
only integrals of the exact solution, but they are also exactly con- 
served by the leapfrog mapping. On the other hand, the extrapola- 
tion procedure does not necessarily conserve this relation. 

The transformation from the variables X to centre-of-mass 
coordinates r can be done as 



ri = 



r-2 



r 3 



(m 3 X 2 — m 2 X 3 ) 
M 

(m 1 Xz — m 3 Ii) 
M 

(7712X1 — 777.1X2) 
M 



(9) 



and the velocities obey the same rule. 



4 INITIAL CONDITIONS 

The systems we studied were equal mass systems. We set all the 
masses m k — 1 (also the gravitational constant was chosen to be 
G= 1). 

We attempted to produce random initial conditions although it 
is not clear how it should be done since the phase space is infinite. 
Initial experiments with various ways showed that the results con- 
sidered in this paper do not depend much on how the initial values 
are randomised. We finally chose the following algorithm: 

First we selected all the 9 coordinates x k from a uniform dis- 
tribution inside a 9 dimensional sphere and similarly the velocity 
components v k . Those were reduced to centre-of-mass system and 
scaling factors c r and c v for the coordinates and velocities respec- 
tively were determined from the conditions 



E = -1, and T/U = k, 



(10) 



where T, U and E — T — U are the kinetic energy, potential and to- 
tal energy respectively, and k, the virial ratio, is a selected quantity. 
After obtaining the scaling coefficients we simply replace 



(ID 



We considered the three values k = 0, k = 0.1 and k = 0.5. The 
first one thus gave a system starting at rest (often called 'the free- 
fall problem') and thus being a plane problem with zero angular 
momentum. The second value k = 0.1 gives systems with small 
initial velocities and the value k = 0.5 corresponds to systems 
starting at the state of virial equilibrium although here it only means 
larger initial velocities and non-negligible angular momentum. 



5 RESULTS 
5.1 Experiments 

We conducted 10,000 experiments for each value of the initial virial 
coefficient k — 0, 0.1, 0.5. Thus altogether 30,000 simulations 
were run to the point when the system disrupted, or the system 
time exceeded 100,000 time units. There were only 383 cases in 
which the computation was halted because of exceeding the time 
limit. There was no reason to halt any of the experiments due to 
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Figure 2. The scatter diagram of disruption times (td) and Lyapunov times 
(te) for the free-fall experiments. The line is a median line determined such 
that there are equally many dots under and above the curve at every value 
of t e . Some systematic initial value effects are visible at small times. 



Figure 3. The scatter diagram of disruption times (td) and Lyapunov times 
(t e ) with the median line. 



loss of accuracy. Our measure of accuracy, the ratio of energy error 
SE and the Lagrangian L — T + U satisfied 



\SE\/L < 1.36 • 10" 



(12) 



in all cases (the number quoted is the maximum), while most values 
(in 24350 cases) were less than 10~ 13 . 





5.2 Scatter diagrams 

Figure |2]illustrates the scatter diagram of (t e ,td) (times according 
to Eq.(|2j) in logarithmic scale for the free-fall experiments. The 
line shown is a median produced by sorting the results according 
to t e , then dividing the whole range into 50 intervals each contain- 
ing equally many cases for which the median was determined by 
sorting this group according to the disruption times td- 

Similarly, Figure [3] illustrates the (t e ,td) scatter diagram of 
all the other experiments. The similarity is obvious although the 
scatter is somewhat larger. 

A third way of illustrating the correlation of (t e ,td) is given 
in Figure [4] where the medians are plotted in the same figure. Since 
in Figures[2]and[3]most cases are concentrated in a much narrower 
region than the total scale, we restricted this figure into the interval 
.94 < t e < 35.2. These limits were determined such that 5% 
of the cases are in the region t e < .94 and similarly 5% in the 
region t e > 35.2. Thus a majority (90%) of data are included in the 
shown interval. Another reason for this restriction is that, at small 
times, the systems disrupt almost immediately and the results are 
likely to be affected by initial value selection. The largest times are 
affected by long ejections without escape, an effect that may also 
partly result from initial value sampling. 

Here three different median curves are produced: one for the 
free-fall experiments, one for small and medium angular momen- 
tum and one for large angular momentum. These were produced by 
sorting the experiments according to the angular momentum value, 
the smaller half of these were included into the small and medium 
group and the rest were considered to have large angular momen- 
tum. One can see that the large angular momentum curve is some- 
what above the others, but the difference is minor. The straight line 
in the figure represents the relation td = 10tJ' 8 , which roughly fits 
the general trend in this region which contains 90% of the data. On 
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Figure 4. Median lines for zero, small and large angular momentum sys- 
tems. The straight line is for t d = 10 t\^ and the circles for t e = 
min(5 tg' 3 , 62 t e ). The curve above others is for the large angular momen- 
tum experiments. The two others, for zero and small angular momentum, do 
not differ noticeably. 

the other hand, the median curves could be better fitted with two 
lines: one with td ~ 5tjj' 3 for t e < 6 and with td ~ 62t e for 
t e > 6, i.e. for large times td oc t e (circles in the figure). 

5.3 Marginal distributions 

It may be of some interest to see what are the marginal distributions 
of the various quantities studied here. We found that, for large val- 
ues of the times, both td and Z = In \5r | = td/t e have exponential 
probability densities which are approximately given by 

ip(t d ) « aexp(-crt d ), a = 1/250 (13) 
tp(Z) « /3exp(-pZ), j3 = 1/45. (14) 

These results are illustrated in the Figures [5] and [6] If one assumes 
those approximations to be valid for small values also then it is easy 
to derive the probability density for t e as 

ip(t e ) = J S(t e — td/Z)aexp(—atd) exp(—/3Z)dt e dZ 
= af3/(ate + P) 2 , (15) 
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Figure 5. Probability density of the disruption times. The line represents 
the function ctexp(— at^) with a = 1/250 
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Figure 6. Probability density of Z = td/te i.e. the disruption time in 
units of the Lyapunov time. The line drawn represents /3exp(— PZ) with 
P = 1/45. 
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which is qualitatively right, as can be seen from Figure [7] where 
the probability density of t e is illustrated. Numerically this result is 
not quite correct, presumably because the exponential functions for 
the probability densities of td and Z are only valid for large values. 



6 CONCLUSION 

The main result in this investigation is the obvious existence of a re- 
lation of the Lecar-Franklin-Murison type in the general equal mass 
three-body problem. In the region where most (90%) of the cases 
in our diagram are located the median disruption time can be rather 
well approximated by a p ower law of exponent n sa 1.8 which 
happens to be close to what lLecar. Franklin. & Murisonl dl992l) ob- 
tained in their pioneering paper. This is somewhat surprising since 
the systems are totally different. On the other hand, our results 
could also be interpreted in the way that for t e < 6 the exponent is 
n ~ 2.3 and for larger t e the exponent is simply n ~ 1. If this is the 
case, one can conclude that for large t e , i.e. for weakly chaotic sys- 
tems, the disruption time is typically proportional to the Lyapunov 
time; 

iLecar et a i] d200lh suggested that the relation in Eq.lQ} is valid 
in Solar System in regions where there are overlapping mean mo- 
tion resonances. However, in the equal mass three-body problem 
there are hardly any resonances causing the behaviour, thus the phe- 
nom enon may be more direct l y relat ed to chaos than to resonances. 
Like lMorbidelli & Froeschiel Jl996l) . we too do not have any theo- 
retical explanation for this phenomenon, which seems to be more 
common than expected. 

In addition to the above, it is interesting to note that (for large 
times) the system lifetimes td have exponential probability density 
both when measured directly in terms of the systems dynamical 
time and also when measured in units of the Lyapunov time t e . 
A consequence of this is that the Lyapunov times have a different 
probability density ipit e ) cx t~ 2 . It would be most interesting to 
find a theoretical explanation for these facts, but we are unable to 
provide one. 
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Figure 7. Probability density of the Lyapunov times t e . This density is 
clearly not exponential but, for large t e , is close to the function 1.8/tg 
drawn in the figure. 



